The Stochastic Green Function (SGF) algorithm 
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We present the Stochastic Green Function (SGF) algorithm designed for bosons on lattices. This 
new quantum Monte Carlo algorithm is independent of the dimension of the system, works in 
continuous imaginary time, and is exact (no error beyond statistical errors). Hamiltonians with 
several species of bosons (and one-dimensional Bose-Fermi Hamiltonians) can be easily simulated. 
Some important features of the algorithm are that it works in the canonical ensemble and gives 
access to n-body Green functions. 
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I. INTRODUCTION 



In the last twenty years, numerical methods have 
gained importance with the increasing power of comput- 
ers. They have been solicited in situations where ana- 
lytical results are missing, due to the complexity of the 
studied problems, and where approximate methods fail to 
give a correct description. But even with computers, ex- 
act calculations are often limited to special cases. This is 
especially true for quantum many-body problems where 
the size of the Hilbert space grows exponentially with 
the size of the system, restricting exact diagonalizations 
to very small systems. Quantum Monte Carlo (QMC) 
methods have been developed in order to simulate bigger 
systems, and allow a correct description of quantum fluc- 
tuations, which are usually missed by mean-field theories 
[l|. QMC methods have given rise to various kinds of 
algorithms @, i 1 i, i, 0,1, i, 0, [Hi • 

We propose here a new algorithm designed for bosons 
on lattices, the Stochastic Green Function (SGF) algo- 
rithm. The algorithm is independent of the dimension of 
the system, and works in continuous imaginary time. The 
algorithm is exact, in the sense that it has no error be- 
yond statistical errors. Hamiltonians with several species 
of bosons Il2. lip. Il4ll (and one-dimensional Bose-Fermi 
mixtures [la. Ila. Il7l|) are easily treated. An important 
property is that the SGF algorithm works in the canoni- 
cal ensemble. This is especially important when working 
with several species of particles, because it is numeri- 
cally difficult to control several numbers of particles in 
the grand-canonical ensemble. Indeed, working with sev- 
eral species in the grand-canonical ensemble requires one 
chemical potential per species. Those chemical potentials 
have to be tuned to flnd the desired number of particles. 
But the difficulty comes from the fact that the number 
of particles of a given species depends on all chemical po- 
tentials. Working in the canonical ensemble makes things 
much simpler, by just choosing the number of particles 
for each species. Another property of the algorithm is 
that it provides access to n-body Green functions, allow- 
ing the calculation of momentum distribution functions 
which are important for the connection between theory 
and experiments. 



II. MOTIVATIONS 



We consider a Hamiltonian of the form 



n = v-T, 



(1) 



where V and T are respectively diagonal and non- 
diagonal (positive definite) operators. We would like 
to have at our disposal an algorithm that is simple and 
able to simulate any Hamiltonian of the form ([T]), in the 
canonical ensemble (for the reason above). As an exam- 
ple, we will consider a Hamiltonian describing two species 
of particles, atoms and diatomic molecules. The parti- 
cles from one species interact together, and interact with 
particles of the other species. We also take into account 
the possibility for atoms to be converted into molecules, 
and vice-versa. The situation can be described by the 
following V and T operators: 



V = 



r 



Uaa Y. "^° i^^ " l) + f^"™ Y. ' 
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1) 



(2) 



gY (™Ia^a» + aUl^n,) (3) 



V and T correspond respectively to the potential and ki- 
netic -|-conversion energies. The aj and a^ operators {mj 
and TOj ) are the creation and annihilation operators of 
atoms (molecules) on site i, and nf — ala,j^ [h™- — mlm^) 
counts the number of atoms (molecules) on site i. The 
sum (j,i) is over pairs of first nearest neighbors. We can 
see that the T operator allows atoms and molecules to 
jump onto neighboring sites, and that two atoms can be 
transformed into one molecule (and vice- versa). The to- 
tal number of atoms Na = X^i ^t ^^"^ the total number of 
molecules N^ ~ X^i ^T ^'^s ^'^^ conserved, however the 
number N = Na + 2iVm is conserved. This is our canon- 
ical constraint. We will not discuss at all the physics of 
this Hamiltonian, referring the interested reader to the 
literature [il,[il[il. 



Among existin g Q MC algorithms, the Canonical Worm 
(CW) algorithm [I^, [li| is a good choice if one wants to 
work in the canonical ensemble and to have access to 
Green functions. However this algorithm makes use of 
a "worm operator" VV, and some complexity of the al- 
gorithm arises when T does not commute with W (see 
section III.B). It is always possible to make the trivial 
choice yV = 1 + T for the worm operator which leads 
to a zero commutator, [VV,'/'] = 0. But such a choice 

is not appropriate at all when the T operator connects 
only neighboring sites (which is usually the case). In- 
deed this would lead to a worm operator that is unable 
to generate spatial discontinuities of the worldlines for 
which the broken parts are separated by more than one 
lattice site. Therefore it would be impossible to measure 
Green functions which require long range discontinuities 
of the worldlines. Moreover, this choice for the worm 
operator generates only local updates, which are known 
to be much less efficient than global updates. Finally, 
those local updates cannot sample the winding, which is 
a quantity of interest when working with periodic bound- 
ary conditions. As a result, a more complicated choice 
has to be made for W. For our chosen T operator, it is 
not trivial to find a suitable W operator that commutes 
with T and satisfies the requirements just mentioned. 
While a suitable W operator might exist, we have not 
managed to find one that can be easily handled. Ref- 
erence [10] proposes an extension of the applicability of 
the worm operator, but this goes beyond our purposes 
of simplicity and generality. The SGF algorithm we pro- 
pose is an alternative way to simulate any Hamiltonian 
of the form (1) in a very simple and general way. Basi- 
cally, once one has a SGF computer code that simulates 
a given Hamiltonian, the only thing to do to extend the 
code to another Hamiltonian is to change the definition 
of the Hamiltonian in the code. 



III. THE ALGORITHM 

A. The partition function and the "Green 
operator" 

The SGF algorithm is derived from the CW algo- 
rithm. We start by considering the partition function 
Z{P) = Tr e~^'^ , and we perform the expansion 



Z{I3) - Tre-^'^T.e/o'^^")"^ 



-l-oo „ 

= Tr ^ /e-'^'^t(T„) • • • t(T2)r(Ti)dri • • • dT„,(4) 

where T,- is the "time ordering" operator and T{t) is 
defined by 



By introducing complete sets of states / — ^^ |V')(V'| 
between each non-diagonal operator T, we get 

Z{P)=Y, feo|e-'^^r(T„)|V'n-l)(^„-l|r(T„_i)|V'„_2) 

X •••(^fc|r(Tfc)|V'fc-i)--- (6) 

X (V'2|t(T2)|V'l)(V'l|^(n)|V'0)dTl---dT„. 

Using the notation Vk for the eigenvalue of V in the eigen- 
state IV'fc), Vk = (^fcjVlV'fc), each matrix element in ^ 
takes the form 



(^,|r(r)|V/) 



^^rV, 



(Vfc|r|^,)e^^^'. 



(7) 



It is useful here to give an interpretation of expression ([6]). 
We assume for the simplicity of this interpretation that 
we have only one species of particles on a one-dimensional 
lattice, and that the T operator is the usual one-body 
operator that makes the particles jump onto neighbor- 
ing sites. The partition function is a sum over all pos- 
sible configurations of time indices ti , • • • , r^ and states 
{IV^fc)}. Figure [T] (left image) shows a representation 
of a possible configuration. We start at imaginary time 
r = with a state h/'o) that contains 3 particles. Then 
the state evolves with the operator e~'^^^° until time ti. 
During this evolution, the state does not change because 
the V operator is diagonal. At time ri a T operator acts 
onto the state, leading to a sum of several new states. 
In this sum of states, only the state jf/'i) survives when 
making the scalar product with the bra (-01 1 . This new 
state differs from k/jo) by a jump of only one particle, 
since we have assumed in our example that T is a one- 
body operator. Thus at time ri one particle jumps onto 
a neighboring site. The new state Wi) then evolves with- 
out changing with the operator e^t'^^^'^i'^i until time T2. 
At time T2 one particle jumps onto a neighboring site 
leading to the new state V'2)--- and so on, until time r„ 
where a last jump of one particle leads to the initial state 
h/io), which evolves without changing with the operator 
g-(/3-r„)y„ m^til time (3. As a result, one configuration 
of time indices ri, • • • , r„ and states { hAfc)} corresponds 
to a set of lines (the worldlines) that the particles follow. 
Because the partition function is a trace, the same state 
appears both at the begining and the end of the imag- 
inary time evolution: The worldlines are periodic with 
period j3. So the partition function has been written as 
a path integral. 

In order to sample the partition function ©, we de- 
fine an extended partition function Z{/3,t) by breaking 
up the propagator e~^^ at imaginary time r and intro- 
ducing a "Green operator" Q, 



Z{P,t) 



Tr e 



-(/3- 



'^"^Ge 



-tH 



(8) 



It is straightforward from ^ to show that the extended 
partition function Z{f3,T) takes the form 



tV^„-tV 



T{t) = e^'^Te 



(5) 
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FIG. 1: Representation of a given configuration of time in- 
dices Ti, ■ ■ ■ ,T„ and states { h/'fc ) I of I'he partition function 
([6]) (left image) and the extended partition function ([9| (right 
image) . 

X ■ ■ ■ {^L+l\fiTL)\i^L){i'L\GiT)\i;R){^R\fiTR)\^R^l)i9) 
X •• ■{lp2\T'{T2)\lpl){lJjl\'t{Ti)\lpQ)dTi ■ ■ -dTn 

where we denote by \4'l) and r^ (h/'i?) and tr) the state 
and the time of action of the T operator appearing to the 
left (right) of the Green operator, and G{t) is defined by 



Gir) 



.tV 



'Ge 



-tV 



(10) 



In order to define the Green operator G , we first introduce 
the " normahzed" creation and annihilation operators A^ 

and A^ 



it = 



1 



Vn + T 



A = 



1 



VnTT 



(11) 



where a' and a are the usual boson creation and anni- 
hilation operators, and n = a^ a is the number operator. 
While unusual, the number operator fi appearing in the 
denominator of a square root is perfectly well defined by 
a power series, 



1 



E 

p=0 



^(2p-l)!! 



p\ 



y/fi + l 
It follows from HIl) and ^^ that 

A^\n) = \n + l) A\n) = \n-l), 



(12) 



(13) 



with the particular case that -4|0) = 0. Apart from this 
exception, the operators A^ and A change a state |n) by 
respectively creating and annihilating one particle, but 
they do not change the norm of the state. 

Using the notation {ip\iq } to denote two subsets of site 
indices ii,i2, • • • ,?p and ji,J2, ■ ■ ■ ,jq with the constraint 
that all indices in subset i are different from the indices in 
subset j (but several indices in one subset may be equal) , 
we define the Green operator G by 



^=ee5p. e n^.n^' (14) 



where gpq is a matrix that will be defined later (see sec- 
tion III.C.6). The Green operator can be viewed as a gen- 
eralization of the "worm operator" introduced in the CW 
algorithm (see section III.B). Note that, because the two 
subsets {ip} and {jq} have no index in common, there 

is no possible cancellation between the operators .4^ and 
A appearing in p4|) . This Green operator is going to be 
sampled stochastically, each configuration leading to a 
measurement of a randomly selected n-body Green func- 
tion, thus justifying the name of the algorithm. 

Let us now consider a state lipL) which is obtained 
from a state \4'r) by creating p particles on sites {ipj 
and destroying q particles on sites {jq}- From (fT4|) we 
can get the corresponding matrix element of G, 



{lpL\G\tl^R) = gpq- 



(15) 



1=1 



In particular, all diagonal matrix elements ("(A^hA) are 
equal to goo, which we will set to unity. The interpre- 
tation of the extended partition function Z{/3,t) is the 
same than the partition function Z{j3), with the addi- 
tion at time r of the Green operator. In the example 
of figure [T] (right image), the Green operator makes two 
particles jump. Note that these jumps are not restricted 
to neighboring sites, in our example one particle jumps 
onto a neighboring site and the other jumps onto a second 
neighboring site. 



B. The update scheme 

As in the CW algorithm in which the worm operator 
updates the configurations of the partition function, we 
use the Green operator to update the configurations ap- 
pearing in 0. But the procedure we follow is different 
and simpler. More precisely, in the CW algorithm the 
worm operator W{t) suggests to create a new T opera- 
tor at time r. This creation is always possible. Then a 
time shift At of the worm operator is chosen, to the left 
or to the right. If the worm operator meets a T oper- 
ator, then it tries to destroy it. This destruction is not 
always possible. When it is not, then the worm operator 
tries to "pass" the T operator. After succeeding to pass 
the operator, a new time shift is chosen and the worm 
keeps moving until reaching another T operator, or until 
the chosen time shift is exhausted. The "passing" pro- 
cedure is always possible only if the commutator of the 
worm operator and the T operator is zero, [yV,T] — 0. 
If it is not, then it will sometimes occur that the worm 
operator cannot pass, and the update will have to be 
cancelled, which leads to some complexity of the algo- 
rithm. In particular all changes made in the operator 
string from the begining of the move must be recorded 
in the event of the need of a restoration. Moreover, it is 
no longer guaranteed that the algorithm is ergodic. In- 
deed, when a rejection occurs because of the unability 
to pass an operator, this rejection is systematic (the 
move is always rejected for the considered configuration) 



instead of statistic (it has a probability to be accepted 
or rejected). This might cause problems with ergodicity. 

In the SGF algorithm, this difficulty is overcome 
thanks to the Green operator. The definition of the 
Green operator ensures that it is always possible to de- 
stroy a T operator. As a result neither a "passing" proce- 
dure nor a zero commutator between Q and T is required. 
In this way the algorithm is simpler. 

The update scheme is the following: 

• We choose a direction of propagation "left" or 
"right" for the Green operator, according to some 
probabilities P{^ and P{~f)- 

• We chose with a probability -PJ_(t) (or Pj^(r)) to 
create a new T operator at time r on the right (or 
the left) of the Green operator. 

• If the creation is accepted, a new intermediate state 
|V') is chosen with some probability P{ip). 

• Then we choose a time shift Ar with a probabil- 
ity P^{At) or P_^(At). If the time shift can be 
exhausted without reaching a T operator, then the 
Green operator is shifted to the new position and 
the update stops there. 

• If a T operator is met before the end of the shift, 
it is destroyed and the Green operator stops there. 

By creating and destroying T operators, the time in- 
dices Tk and states lipk) visited by the Green operator 
are updated. This way the extended partition function 
2{P,t) ^ is sampled. When a diagonal configuration 
of the Green operator occurs, KAl) = KA-r)) such a par- 
ticular configuration of Z{(3,t) belongs to the space of 
configurations of Z{l3). Measurements of physical quan- 
tities can then be performed. Since it is always possible 
to create an operator at any time, and since it is always 
possible to destroy a reached operator, it follows that the 
ergodicity of the algorithm is ensured. 



C. Detailed balance 

We describe here how to perform the update scheme 
by satisfying detailed balance. Four different situations 
have to be considered (versus five in the CW algorithm, 
the extra one being the "passing" move). 

1. No creation, shift, no destruction. 

2. Creation, shift, no destruction. 

3. No creation, shift, destruction. 

4. Creation, shift, destruction. 

We will assume in the following that a left move is chosen. 
We denote the probability of the initial (final) configura- 
tion by Pi {Pf). We call Si^f the probabihty to suggest 



a transition from configuration i to configuration /, and 
Sf-+i the probability of the reverse transition. Finally we 
call Ai^f the acceptance rate of a transition from i to /, 
and Af^i the acceptance rate of the reverse transition. 
The detailed balance can be written 



P,S,^fA^f 



PfSf^^Af^,. 



(16) 



A possible solution for the acceptance rate is the 
Metropolis solution [ISj . 



ith 



Ai^f = min ( l,g 



PsSl-,^ 



(17) 



(18) 



1. No creation, shift, no destruction 

We consider here the case where a left move is chosen 
with probability P(<— ), no creation is performed with 
probability 1 — Pjl (r) , a time shift of Ar is chosen with 
probability P^(Ar), and finally no destruction occurs 
because the time shift is supposed to be too small to 
reach a T operator. 

The probability of the initial configuration is the Boltz- 
niann weight appearing in ^: 



P^ (X {tpL\Q{r)\'rPR) 
The probability of the final configuration is: 



(19) 



Pf oc (V'L|^(T + Ar)|Vfl) 

ex e(^+^^)^-(VL|^|^fl)e-("+^"'^« (20) 

The probability to suggest the transition from the initial 
configuration to the final configuration is the probability 
P(<— ) to choose a left move, times the probability of no 
creation 1 — PJL(t), times the probability P^(Ar) to 
perform a left shift of Ar: 



s.^f = p{^){l-pl{T))p^{^T) 



(21) 



The probability to suggest the reverse move is exactly 
symmetric: 



Sf^, = P(^)(l -PI{t + Ar))P^iAr) 



(22) 



The acceptance rate of the corresponding move is given 
by (Ull), with 



_ e-^-(^^^^«)P(^)(l - PUr + Ar))P^(Ar) 
'^ P(^)(l_pi(r))P^(Ar) 

Because of the exponential appearing in (P5)l the accep- 
tance rate might be small if the diagonal energy Vr is 



greater than Vl- In order to keep a good acceptance 
rate, this exponential can be cancelled by making a good 
choice for the probability of the time shift, 

-ArVn p_^(Ar) = Vie-^^^^ (24) 



P^(At) = Vrc- 



Equation 



becomes 
Vl 






^'' P{^)il - Pl{r)) 

where we have defined t' — t + At and V^ = Vr, and 
we have used the notation q^^ to emphasize that there is 
no creation and no destruction. We also have explicitly 
writen q^^ as a product of a quantity that depends only on 
the initial configuration, times a quantity that depends 
only on the final configuration. 

2. Creation, shift, no destruction 

We consider here the case where a left move is chosen 
with probability P(<— ), a creation of a T operator is per- 
formed with probability PJL(t) thus introducing a new 
intermediate state tpjir \ on the right of the Green oper- 
ator, chosen with a probability P{ipiii). A time shift of 
At is then chosen with probability P^(Ar), and finally 
no destruction occurs because the time shift is supposed 
to be too small to reach a T operator. 

The probability of the initial configuration is given by 
(fT9)) . The probability of the final configuration is: 



Pf(x{^L\giT + AT)\^R,){i;R,\fiT)\^R) 

The probability to suggest the transition from the initial 
configuration to the final configuration is the probability 
P(^) to choose a left move, times the probability PJL(t) 
to create a new T operator at time r, times the prob- 
ability P{tpRi) to choose the new state IV'fl/): times the 
probability P^{At) to perform a left shift of At: 

S^^f - P{^)pI{t)P{4,r.)P^{At) 



= P{^)Pl{T)P{^R,)Vi 



,-AtV„ 



R,e -■•- (27) 

We have explicitly used our previous choice ([M]) for the 
probability of the time shift, by taking care that the state 
Wr) on the right of the Green operator has been updated 

to l^fl/). 

The probability to suggest the reverse move is the 
probability to choose a right move P{-^), times the prob- 
ability of no creation 1 — Pj^ (r -I- At) , times the probabil- 
ity to reach the T(r) operator on the right and destroy 
it. This latter probability is the probability to choose a 
right shift greater than At. It is obtained by integrating 
P_^(t) from At to -|-cx). Because of our choice (PU, this 
integral can be explicitly calculated: 

p + CC 

Sf^., = P(^)(l-Pl(r + Ar)) I P^{t)dt 



At 



P(^)(l-Pl(r'))e- 



Using (|18p to calculate the corresponding acceptance fac- 
tor qc^, all exponentials cancel and we get: 

_ (V^L|g|V'flO(V'fl>|r|v>H)P(^)(i-Pl.(rO) 

'^'^ {i^L\G\i^R)P{^)PUr)P{^R,)VR, 

We can here explicitly make a choice for the probability 
P{'ipR') of the new state \ipR')- If we choose the new 
state proportionally to the Boltzmann weight of the new 
configuration. 



Pi^R') 



{'4'l\G\^r'){iPr'\T\iPr) 

S^H' {'^L\G\lpR'){lpR'\f\lpR) 
{i>L\G\ll'R'){lpR'\'t\ll^R) 



{i^L\GT\^PR) 

then the acceptance factor (|29p becomes 

{ijL\Gf\tljR) p(^)(l_pl^(r')) 



(30) 



qcfi = 



P(^)Pt(T)(V'L 1^1 Vfl.) 



Vr, 



(31) 

where q^^ is written as a quantity that depends only on 
the initial configuration, times a quantity that depends 
only on the final configuration. 



3. No creation, shift, destruction 

We consider here the case where a left move is chosen 
with probability P(<— ), and no creation is performed with 
probability 1 — Pjl (r) . A time shift of At is then chosen 
with probabihty P^(At), and a destruction of the T{tl) 
operator to the left of the Green operator occurs because 
the chosen time shift At is taken to be larger than t^ — t. 

The probability of the initial configuration is: 



Pi(x{^pL+l\f{TL)\^pL){■^^^L\G{T)\^^^R) 



cxe 



tlVl^ 



'TVr, 



(32) 



The probability of the final configuration is: 

Pf oc (VL + l|^(TL)|Vfi) 

ex e"^^^+i(V;L+i|^|7/.fl)e-^^^« (33) 



-AtVl 



(28) 



The probability to suggest the transition from the initial 
configuration to the final configuration is the probability 
P(^-) to suggest a left move, times the probability of no 
creation 1 — P,L(t), times the probability to suggest a 
shift greater than tl — t, that is to say the integral of 
P^{t) from Ti — T to +oo: 

S.^f - P(^)(1-P1(t)) / P^it)dt 

= P(^)(l-Pl(T))e-(^^-^)^« (34) 

The probability of the reverse transition is the probabil- 
ity to choose a right move, times the probability P1^{tl) 



to create a new T operator to the left of the Green oper- 
ator, times the probabUity P{'ijj]^) to choose the new m- 
termediate state |V'l)i times the probabihty P^{tl — t) 
to perform a right shift of tl — t. We get 



Sf^^ = P{^)pI{tl)P{^/jl)P^{tl ~ t) 



with 



P{^L 



{tljL + l\T\lpL){lpL\G\lpR) 



(35) 



(36) 



{ipL+i\Tg\^R) 

The acceptance factor is given by 

Vl ^, P{-^)Pl{T'){^PL'\g\^R') 

^''' Pi^){l-Plir))'' {iJL'\fg\i^R') 

(37) 
where we have defined r' = r^, lipL') — hAi+i/i f-nd 

V'-R') — ■0-r)- Again, the acceptance factor q^d is writ- 
ten as a quantity that depends only on the initial config- 
uration, times a quantity that depends only on the final 
configuration. 

4- Creation, shift, destruction 

Finally we consider the case where a left move is chosen 
with probability P(<— ), a creation of a new T operator is 
chosen with probability PjL(r), leading to the introduc- 
tion of a new state h/'-R') with probability P('0_r')j ^^^ 
the T operator to the left of the Green operator is de- 
stroyed because the chosen time shift At is taken to be 
greater than tl — t. 

The probability of the initial configuration is given by 
([32)1 . The probability of the final configuration is: 

\G{TL)\lpR'){lljR'\'I'{T)\'ll;ii) 

i(V'L+i|^|V'fl')e-("^-")^«'(Vfl'|t|Vfl)e-"^'^(38) 

The probability to suggest the transition from the initial 
configuration to the final configuration is the probability 
P(^) to chose a left move, times the probability PJ1(t) 
to create a new T operator to the right of the Green 
operator, times the probability P{'ipf{i) to choose the new 
state IV'fl'): times the probability to reach the T to the 
left of the Green operator, that is to say the integral of 
P^{t) from Ti — T to -l-oo: 

r + co 

s^^f - p(^)Pl(T)P(VflO / p^it)dt 

= P{^)Pl{T)P{ijR,)e~^^^-^^^'^' (39) 

The probability to suggest the reverse transition is ex- 
actly symmetric: 



P/(x(V'L + l 



cue 



Sf^., = P{-^)PI{tl)P{'^l) / P^{t)dt 

= P(^)Pl(ri)P(VL)e-(^^-^)^^ (40) 



Using the notation r' = tl and \4'l') = |V'l+i), the 
acceptance factor takes the form 



qcd = 



{iPl\GT\iPr) 



P{^)PI{t){^l\G\^r) 



P{-^)Pl{r'){i>L'\G\i^R') 

{lpL'\'tG\lpR/) 

(41) 

and is again written as a quantity that depends only on 
the initial configuration, times a quantity that depends 
only on the final configuration. 



5. Simplification of the acceptance factors 

Having determined all acceptance factors 
(lgf<f,qc4,q,/d-,qcd for aU kinds of updates, we still 
have some freedom for the choice of the probabilities 
of creation P^{t) and P1^{t), and the probabilities of 
choosing a left or right move, P(<— ) and P(^). 

Let us define the following quantities: 



9^(r) = 

q^ir) = 



{i^L\Gf\i^R) 



P{^)Pl{T){ijL\G\i^R) 

{'>Pl\T'G\-iPr) 

P(^)PX(r)(VL|^|Vfl) 

Vl 
P{^){l-Pl{r)) 

Vr 
P(-)(1-Pi(r)) 

The acceptance factors take the form: 



q</d = 



qAr') 



qc^' 



qcd 



qMr') 



(42) 

(43) 
(44) 
(45) 

(46) 

(47) 



We immediately see that all acceptance factors become 
equal if q^ir) — 4^{t) and q^ir) = i^ij). This is 
realized if we choose for the probabilities of creation: 



Pli.r) 
Plir) 



{iPl\GT\iPr) 



VL{^L\G\i}R) + {^L\GT\i>R) 



VR{^L\G\ijR) + {^Pl\TG\^Pr) 
Then all acceptance factors q^^, qa^, q^d, qcd become 
P{^)r^{T) 



(48) 
(49) 



with 



P{^)r^{T') 
P{-^)r^{r') 



-ir) = VL + 



for a left move (50) 

for a right move, (51) 



{4'l\GT\iPr) 

{iIjl\G\iPr) 



(52) 



r^{T) ^Vr + 



(V'L|^g|V'fl) 



(53) 



Finally, we still have some freedom for the choice of P(<— ) 
and P{-^)- If we choose 



P(H = 



-{r) 



r^ir) + r->{T) 



PH) = 



.W 



r^(T) + ?'_*(t) 



(54) 



and define 



R^ 



r^{T)+r^{T) Rf ^r^{T')+r-,{T'), (55) 



then we are left with a unique acceptance factor which 
is independent of the chosen direction of the move, in- 
dependent of the nature of the move (creation or not, 
destruction or not), and depends only on the initial and 
the final configuration: 



Ri 

R'f' 



(56) 



This result allows to simplify the algorithm. Indeed, by 
combining ([SS]) and P^ we get 



which can be rewritten as 



Ri 
R'f" 



J^iX^iJ^i 



= 1. 



*/ 



(57) 



(58) 



This last equation can be interpreted as follows: If we 
accept all moves without taking care of the acceptance 
factor, then we are sampling the extended partition func- 
tion according to the pseudo Boltzmann weight Pg = RP. 
The algorithm is simplified, because all moves are ac- 
cepted and it is not necessary to keep track of all changes 
performed during an update in case of the need of a 
restoration of the initial configuration. 

The statistics of a physical quantity described by an 
operator O relevant to the real Boltzmann distribution 
is recovered by using the relation 



(O), 



(o/R) 



-Ps 



(i/i?), 



(59) 



which is well defined because the quantity R is well be- 
haved: it never vanishes nor diverges. We emphasize here 
that this simplification of accepting all moves is always 
possible in the SGF algorithm, even if the T operator 
does not commute with the Green operator, whereas it 
is possible in the CW algorithm only if the T operator 
commutes with the worm operator. 



6. Determination of the g^q matrix 

Measurements of physical quantities represented by di- 
agonal operators can be performed only when the Green 



operator is in a diagonal configuration, wl) = Wr)- 
The situation is different when measuring Green func- 
tions: Their measurement extends from a diagonal con- 
figuration to another, while exploring the extended space 
of configurations. But the end of the measurement is 
marked by the return back to a diagonal configuration 
(see section III.D). For the Green operator to have a 
chance to go back to a diagonal configuration, an ap- 
propriate choice of the gpq matrix must be done. As a 
result, Qpq must decrease sufficiently fast as p and q go 
to infinity, in order to prevent the state V'l) to be too 
different from \iI)r)- The exact choice of gpq depends on 
the application of the algorithm. It is natural to choose 
gpq to be a decreasing function oi p + q. 

For the example of the Hamiltonian described by ([T]) , 
(m, and ((21), we find that the choice 



9pq 



1 



,-i{2-p-qf 



if p 
if p 



q<2 
q>2 



(60) 



leads to a very good statistics for one-body Green 
functions of the form (a\aj\ or (rnlmj). However if 
one is interested in more complicated Green functions, 
(a\a,akai) for instance, the choice 



9pq ~ 



-i{A-p-qf 



if p 
if p 



g < 4 
g > 4 



(61) 



is more appropriate, accompanied by a slowing down of 
the algorithm but also by an improvement of the statis- 
tics. An important thing to notice is that there cannot 
be any " cutoff' on gpq\ Indeed, the choice 



9pq 



if p 
if p 



g < 4 
g > 4 



(62) 



leads to a crash of the algorithm, because configurations 
where Wl) and Wr) are connected by p -I- g = 4 cre- 
ations and annihilations will occur, and the Green opera- 
tor might be unable to destroy an operator, if its destruc- 
tion leads to states ji/'i) and \'^fj) that are connected by 
a higher order of creations and annihilations (see section 
IV for a concrete example). 



7. Efficiency and purposes of the algorithm 

The generality of the SGF algorithm could result in 
loss of efficiency compared to other algorithms (when 
such algorithms can be applied) because the extended 
space of configurations which is sampled is much larger 
than the one sampled by other methods, due to the in- 
finite sum in the expression of the Green operator. The 
advantage however is that this gives access to n-body 
Green functions. Any configuration (complicated or not) 
of the Green operator allows a measurement of the cor- 
responding Green function (see section III.D). Hence the 
purpose of the SGF algorithm is not to compete with the 
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speed of other algorithms. The properties that make the 
SGF method useful are: (i) it is simple to apply to any 
Hamiltonian of the form (1), (ii) it is very general and 
(iii) n-body Green functions are easily accessible. 



D. Measuring physical quantities 



Let us consider the density operator of the system, 
p — -g-e^'^^. For any physical quantity described by an 
operator O, the expectation value is given by: 



{6) = Trdp 

= ^(Vo|o|V')(V'|p|V'o) 



(63) 



■^oip 



1. Quantities represented by diagonal operators 



If the operator O is diagonal, then (p5)) becomes 



(6) = ^(^|a|v)(V'|/5|^) 

« ^ E '^(^^)' (64) 



tps^P 



where the notation t/js <— p means that the states tps are 
generated according to the Boltzmann weight (^ps\p\'4's), 
and Sd is the number of samples of diagonal configura- 
tions. Equation ()64|) becomes exact when the number 
of samples goes to infinity, and the error decays as the 
square root of Sd, according to the Central Limit theo- 
rem. Since we are actually sampling with a pseudo Boltz- 
mann distribution, equation ([5^]) must be used leading to 



(6) 



E. 



0{i^s)/R{A) 



E^.<-,. yms 



(65) 



where the notation ips <— Ps means that the states ips 
are generated by accepting all moves, irrespective to the 
acceptance factor ([5S)) . As a result, all quantities rep- 
resented by diagonal operators can be directly measured 
when a diagonal configuration of the Green operator oc- 
curs. This includes density-density correlation functions 
(nifij), for instance. In particular, one of the easiest 
quantity to measure is the diagonal energy (V). It is 
measured by averaging the potential Vl (or Vp) to the 
left (or the right) of the Green operator using equation 
([65| . The non-diagonal energy (T) should be evaluated 
in principle by measuring the one-body Green functions 
(described below). But we have actually a direct access, 
simply by averaging the length n of the operator string 



<t> - i(.> 



(66) 



Indeed equation (|66|) can be derived easily by considering 
the quantity Z{P,a) = Tr e-''^*""^). From (g]) this can 
be written as 



Z{P, a) = Tr e 



~PVrj,^^aj^'f{T)dT 



f{T)dT (67) 



Tre-'^'^T,^ 



By noticing that (T) — -hi -^ ln-2(/3, a) ] , we get 



/t) = — ynTre-^*^l 



(3Z 



|q=1 



f{T)dT 



(68) 



Boltzmann weight of n 



which leads to 

We can actually improve the estimates of diagonal 
quantities by integrating them over the imaginary time 
axis. 



{d) = -^j\d{r))dT. 



(69) 



In order to evaluate ([M|) , let us consider a given configu- 
ration of time indices ti , r2 , • • • , r„ of the operator string 
in ([S]), with the convention that t„+i — ti. For any r in 
the range [O, /?[, we have the identity 



^e(Tfe<T<Tfe+l) = l, 



(70) 



fc=l 



with Q{arg) = 1 if arg is true, and otherwise. The 
identity expresses that r has to be located somewhere in- 
between two consecutive time indices Tk and t^+i . There- 
fore we have 



6(7) = a(T)^e(rfc<r<rfc+i) 

k=l 

n 



(71) 



fc=i 



The integral of (|7ip is immediate 



1 ./3 -, " 

- / d{r)dT = -J2 0{rk){Tk+i - Tk) (72) 

The right hand side of ([72|l can be directly averaged over 
the simulation, and leads to an improved estimate of {Oj . 
Time dependent density-density correlation functions, 

C.,(r) - (n,(0)n,(r)) 



= - / n,(T')n,(T + T')dT', (73) 

P Jo 

are also easy to measure using expression ((7T|) . 



The superfluid density ps can be determined by making 
use of Pollock and Ceperley's formula p^ . 



{W')L 



2-d 



2dtl3 



(74) 



where W is the winding number, L is the number of 
lattice sites in one direction of the lattice (assuming the 
same value for all directions), t is the hopping parameter, 
and d the dimension. The winding number is sampled by 
the algorithm, and is easy to measure. It is equal to the 
number of times that the worldlines cross the boundaries 
of the system in a given direction, minus the number of 
times they cross in the opposite direction. This way the 
superfluid density is easily evaluated. Section III.E ex- 
plains how to determine the zero-temperature superfluid 
density, using a finite-temperature simulation. 



2. Quantities represented by non-diagonal operators 

Any physical quantity represented by a non-diagonal 
operator can be expressed in terms of Green functions. 
Green functions can be measured "on-the-fly" while the 
Green operator is updating configurations. Let us con- 
sider the expection value of a particular term Gp of the 
Green operator: 

{Gp) = TrGpP 

= E (V'L|Gp|Vii)(V'«|/5|^L) (75) 

It is important to understand that the states |'0l) and 
\'4'r) are not generated with probability proportional to 
('0_r|/3|'0l) but with probability P('(/;i, i/'fl) proportionnal 
to (V'L|^|V'i?)(V'fi|p|V'L), that is to say 



Pi'^L.i^R.) 



{'4>l\G\iI!r){iI>r\p\'4}l) 



Tr^p 
Thus, equation (|75p can be rewritten as: 

{G,)-T.gpY. %i;#^(^-^-) 



(76) 



i'L,^. 



^ (V'il^lV'fl) 



(77) 



By performing a sampling according to the distribution 
P{i^L,ii'B), we get 



{Gp) 



Tr Qp 
S 

Tr Qp 



■^ {'4'l\Gp\'iPb) 

^ e((V'L|Gp|V'i?)7^0), (78) 



where S is the number of samples including diagonal and 
non-diagonal configurations. In order to evaluate (j78p . 



one needs to be able to calculate Tr Qp. This can be 
achieved by considering the trace of p: 



Tr p 



1 



= E {'4^L\'lpR){i^R\p\lpL) 

1pL,i>R 



■4'L 



-%^ {^L\Q\'^Ii) 



Tr Qp 
S 



^ S{iI'l,iPr) 

IpL^lpR^P 



= f Tr Qp 



(79) 



By injecting (O into ((TS]) we get 

{Gp)^^ E 0((V'L|Gp|^fl)^O). (80) 



Sd 



lf'L,lf'R^P 



Again, since we are sampling by accepting all moves, 
equation ((59|) must be used instead of (|80p leading to 



{Gp) 



E 



l('L;l('R*-Ps 



e((^i,|Gp|v.fl)#o) 
Rii>L,-4>R) 



2.^1/ 



1/R{^jl,^r) 



(81) 



where the notation V'LiV'-R '^ Ps means that the states 
V'l) and h/'-R/ are generated by accepting all moves. Fi- 
nally, a renormalization can be performed onto Gp by 
inverting (jlip in order to get the desired Green function. 
For example, let us suppose that we want to measure 
/a2a5V The corresponding term G25 of the Green oper- 
ator is G25 =511^2^5. We get 

{alas) = {Al^/n2 + l^/n5 + lA5) 
1 



(V^G25V^) 

511 



1 E 



i'L,4'R^P3 



r(l/)I,|A^A5|tAH)A 



RilpL,i'R.) 



511 



Z-^V 



l/i?(7Ai,^fl,) 



(82) 



where 712 and n^ are respectively occupation numbers of 
the states I^Jl) and jV'i?)- 



E. Improved estimator for the zero-temperature 
superfluid density 

As we have seen in the previous section, the superfluid 
density ps can be easily obtained by using (|74p . How- 
ever, the superfluid density shows a strong dependence on 
the inverse temperature /3, especially for one-dimensional 
systems (ID). It has been shown |20| for ID systems that 
superfluidity exists in the thermodynamic limit only at 
zero temperature, and that the zero-temperature limit 
should be taken prior to the thermodynamic limit. This 
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requires in principle to perform simulations with increas- 
ing values of the inverse temperature (3, which is expen- 
sive in computer time, and then perform an extrapolation 
to /? = -|-oo. We propose here an improved estimator that 
gives the zero-temperature superfluid density at arbitrary 
large temperature, thus making simulations easier. 

This improved estimator has been proposed by Ba- 
trouni and Scalettar for a discrete time World Line al- 
gorithm . We give here a generalization to continuous 
time. The improved estimator is actually for the wind- 
ing number, and we determine p^ using ([Ti]) . We con- 
sider here a one-dimensional system, in order to ease the 
explanation of the method. Let us introduce for our pur- 
pose the continuous time pseudo-current j (r) of a given 
configuration of the operator string in ^ , 



n 
j(r)=^D(Tfe)<5(T-Tfe), 



with 



D(rfe) = 



fe=i 



1 if right jump at time r^ 
-1 if left jump at time r^ 



(83) 



(84) 



The winding number is then obtained by integrating the 
pseudo-current over the imaginary time 



W 



j{T)dT 



1 

fc=l 



(85) 



The trick is the following: Instead of directly calculating 
the winding using ()85|) . we evaluate the Fourier transform 
j(w) of dMl) for wi = 27r//3 and wj = 47r//3: 



.?(^) = 



/3 



dr 



/o 

n 

E 
fc=i 



D(rfc)e-'--'= 



(86) 
(87) 



It is straightforward to check that W^ = \j{uj — 0)| /i^. 
But instead of calculating j{lu = 0), we perform an ex- 
trapolation to zero frequency: 



W^ 



2|j('^i)|^- |j(^2)|^ 



/L' 



Wi 



Equation ([88|) becomes exact when /3 goes to infinity, 
since both lji and L02 go to zero. It turns out that, when 
numerically computed, W^^.^^ shows a quasi-linear depen- 
dance in (3. As a result, when injected in (|74|) . the depen- 
dence in (3 is cancelled by the denominator. In this way, 
the zero-temperature superfluid density can be evaluated 
at non-zero temperature. Figure [2] shows the efficiency 
of this method by comparing the dependence in temper- 
ature of the superfluid density, calculated using the true 
winding number W , and using the improved estimator 

Wext. 







' 1 ' 1 
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FIG. 2: (Color online) The superfluid density as a function 
of the inverse temperature /?. Comparison between the value 
ps measured using the true winding number, and the value 
pf^^ measured using the improved estimator. The improved 
estimator converges faster to the large /3 (zero-temperature) 
limit. 



IV. THE ALGORITHM IN PRACTICE 

We describe here how to represent in practice a Hami- 
tonian, the Green operator, and the associated extended 
partition function in the memory of a computer. The 
proposed representation may not be the most efficient, 
but it has the advantage of being easy to handle. We 
will consider here the Hamiltonian ([1]) with the V and T 
operators defined by ^ and ([3]). We have seen that a 
given configuration of the operator string in ^ is fully 
determined by the time indices ri , • • • , r„ and the set of 
states { h/'fc)}- However there is too much information 
in such a representation, because the states in the set 
{li/'fe)} are all almost the same. Thus, it is better to 
specify a configuration by the two states [V'l) and |^_r), 

and specify for each T{Tk) operator which term in ([3]) 
is actually acting. We can use the following "Operator" 
data structure to represent each operator in the operator 
string: 

• Type - An integer number that describes if the op- 



erator is a a] a 



m\m., m\a^a^, 



or alalm- opera- 



tor. A special value is assigned for the Q operator. 

• Time - a real number that represents the time r^ 
of action of the operator. 

• Indexl - An integer number that is the site index: 
If the type is ala,j or mjm,, then Indexl is the index 
of the creation operator. If the type is is m\a^a^ 
or a\a\m^^ then Indexl is the site index where the 
conversion occurs. If the type is Q, then the value 
of Indexl is ignored. 

• Index2 ~ An integer number that is the site index: 
If the type is a\cLj or mjm-, then Indcx2 is the 
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index of the annihilation operator. If the type is 
mla^a^, a\a\m^, or Q, then the value of Index2 is 
ignored. 

• PtrL - A pointer onto an "Operator" data struc- 
ture that represents the operator on the left of this 
operator. 

• PtrR " A pointer onto an "Operator" data struc- 
ture that represents the operator on the right of 
this operator. 

This data structure is part of a doubly linked list. It 
can be used to build the operator string by linking the 
"Operators" together. The states |V'l) and IV'-r) can be 
represented by arrays of occupation numbers. The con- 
figuration of the operator string is then fully represented 
(Fig. [3]) . This structure has the advantage to allow easily 
the insertion of a new piece or the destruction of a piece, 
which correspond respectively to a creation or a destruc- 
tion of a T operator. Changing the time r of the Green 
operator in the range [r^ , t/{] corresponds to moving the 
Green operator between its left and right T operators. 



End 



; 1 G=A,A3M4 rz. v 



T=a|a, 



|V^>=|1011:0000> 
|»|/„>=|0001:0001> 



■^2 



FIG. 3: The operator string can be represented by a dou- 
bly linked list of "Operator" data structures, where each 
piece represents a T operator or the Green operator. The 
advantage of such a representation is that it is easy to in- 
sert or remove T operators. We have used the notation 
|V) = |njn>g< : n'^n'^n^nT) for the states |Vi) and 

It is useful here to add some extra information in the 
computer. We define the "Field operator" data structure 
in order to have a suitable representation of the Green 
operator: 

• Type - An integer describing the type of the nor- 
malized field operator, Aj, A^, Mj, or M^ . 

• Index - An integer describing the site index where 
the field operator is acting. 

• Ptr - A pointer onto a " Field operator" data struc- 
ture that represents the next field operator. 

This data structure is part of a linked list. It can be used 
to build the term of the Green operator that connects 



FIG. 4: The active term of the Green operator can be repre- 
sented by a linked list of "Field operarator" data structures, 
where each piece represents a normalized creation or annihi- 
lation operator. 



the states \'iI^l) and IV'ii:) (Fig- HI)- We will call this term 
the "active term" of Q and denote it by G. 

We have seen in section III.C that we need to be able 
to evaluate matrix elements of the form: 

Nt ^ {ipk+i\f\iPk) (89) 

Ng = {iI'l\G\i1^r) (90) 

NgT = {ll^L\Gf\^l^R) = Y, {ll^L\G\^){'^\f\ll^R) (91) 

iVre = (V-Llt^lV-fl) = ^ (V'L|t|^)(V|^|V'fl,) (92) 



The Nt matrix element is easy to calculate, since we 
know from the " Operator" data structure which term of 
T is acting. The Nq matrix element is also easy to calcu- 
late: We just run over the linked list that represents the 
active term of the Green operator and count how many 
creation operators and how many annihilation operators 
we have. The value of the matrix element is then given 
by the gpq matrix. 

The evaluation of the Nqt (or Ntg) matrix element is 
required when we calculate the probability of insertion of 
a new T operator. For this, we need to look for all pos- 
sible intermediate states k/;). For a given intermediate 
state, only one term of the T operator (for example a\a^ 
or 77120203) gives a non-zero value to the matrix element. 
We will call this term the " active term" of T and denote 
it by T. The important thing to notice is that all active 
terms are inversible and that the inverse of T is propor- 
tional to T^ So the procedure is the following: Instead 
of building the list of states that we get by applying T 
onto \'<Pr) for Ngt (or (V^lJ for Ntg): we build a hst 
of all possible active terms T that give a non zero value 
when applied onto the ket (or the bra). Then for each 
possible active term T we consider the associated nor- 
malized operator T (obtained by replacing all creation 
and annihilation operators by the corresponding normal- 
ized operators (jlip ). and we build the new corresponding 
active term G ' of the Green operator as folows: 

^G" = Grt (93) 

{i>L\Q\'^R) -^ (V'l|7'|V')(7/;|^|V'_r) 

^G'=T'^G (94) 

It is clear that equations ([M)) and ([M]) always have a so- 
lution for G", and that it corresponds to a term of the 
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Green operator. This ensures that it is always possible 
to create a T operator acting on any state at any imag- 
inary time. It may happen that the new active term G ' 
contains normalized creation and annihilation operators 
that cancel each other. In that case a "simplification 
procedure" has to be called in order to remove the obso- 
lete operators and prevent a useless growing of the linked 
list. Having determined all new {G',T) pairs, it is easy 
to calculate the weights of the corresponding matrix el- 
ements. Then a particular pair can be chosen with a 
probability proportional to its weight. A new piece of 
" Operator" data structure is then created and initialized 
with the active term T of the chosen pair, and inserted in 
the doubly linked list of the operator string. The active 
term G of the Green operator is also updated with G '. 

Finally, we need to determine the new active term G ' 
of the Green operator when a T operator is destroyed. It 
is simply given by: 

{'lpL+l\T\tl)L){->JjL\G\-4'R) -> (^L+ll^l-^i?) 

^G'^TG (95) 

{lpL\G\lpR){lpR\T'\lijR-i) -^ (■0L|^|^fl-l) 

^G' = GT (96) 

Again, G ' always has a solution which is a particular 
term of the Green operator. Thus it is always possible 
to destroy any encountered operator. The "simplifica- 
tion procedure" is again called in order to remove ob- 
solete normalized creation and annihilation operators in 
the new G ' . The T operator is removed from the doubly 
linked list of the operator string, and the active term G 
of the Green operator is updated with G ' . 

As a concrete example, let us build the list of all pos- 
sible active terms T of the T operator that can be in- 
serted to the left of the Green operator of Fig. [3l and 
the associated active terms G". We look for all possible 
transitions: 

(1011 : 0000|^|0001 : OOOl) 

-^ (1011 : OOOOJf |V')(V'|^|0001 : OOOl) (97) 
^G' = T'^G 

The solutions (after simplification of G") are: 



a\ a. 



l"'2 



a^a. 



3 "2 



G" 



a^a 



3 "4 



aAQ. 



'4 "3 



AAm, 

AA^i 
a\aIm^ 

aXAAa^^a 
a\a\aIa^m^ 



(99) 



Let us suppose that the (T = a\ai,G ' = A|A|yl|^yl4M4) 
pair is chosen. The new state -0) introduced on the left 
of the Green operator is: 



|0) = A\A\AlAiM4\mQl : Qmi) 
= |2010:0000) 




FIG. 5: (Color online) Comparison between an exact diago- 
nalization on a 4-site lattice, and the SGF algorithm. The 
paramters are ta = 1, tm = 0.5, Uaa = 4, Uam = 12, 
Umm = -|-oo, D = 0, and /3 = 4. The figure shows the total 
energy (S\, the number of atoms (Na), and the number of 

molecules (Nm)- The exact curves fit perfectly in the error 
bars of the QMC results. Note that for all points we have 
Na + 2Nm = 3, which is our canonical constraint. 



Now if we decide to destroy the T operator on the right 
of the Green operator, the state Wr") is removed and 
the only solution for the new active term G" is (after 
simplification): 



G" = G'mIM- 



a\a\aIa4,Mz 



(101) 



The new state \'ip'ji) on the right of the Green operator 
is given by: 



IV'fl) = MIM4\'<Pr) 

= |0001 : 0010) 



(102) 



(100) 



This example illustrates how the algorithm is easy to 
apply to any Hamiltonian of the form ([1]) , provided that 
the non-diagonal part T is positive definite. Figures 
and [6] show a comparison between an exact diagonaliza- 
tion on a 4-site lattice initially loaded with 3 atoms and 
no molecule, and QMC results obtained with the SGF al- 
gorithm. The perfect agreement confirms the exactness 
of the algorithm. 



V. CONCLUSION 

We present a new quantum Monte Carlo algorithm: 
The Stochastic Green Function algorithm. This algo- 
rithm can be easily applied to a wide class of Hamilto- 
nians, including multi-species Hamiltonians. The algo- 
rithm is completely independent of the dimension of the 
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system, and works in the canonical ensemble, which is 
prefered for systems with several species of particles. Fi- 
nally, the algorithm gives access to n-body Green func- 
tions, which provide momentum distribution functions, 
thus allowing useful connections with experiments. 
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